import numpy as np


def kriging(points, values, target_point, model='exponential', nugget=0.0, sill=1.0, range_=1.0):
    n = len(points)

    # 计算距离矩阵
    distance_matrix = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            distance_matrix[i, j] = np.linalg.norm(points[i] - points[j])

    # 计算半变异函数
    def variogram(h, model, nugget, sill, range_):
        if model == 'exponential':
            return nugget + sill * (1 - np.exp(-3*h/range_))
        elif model == 'gaussian':
            return nugget + sill * (1 - np.exp(-3*(h/range_)**2))
        elif model == 'spherical':
            if h <= range_:
                return nugget + sill * (1.5 * (h/range_) - 0.5 * (h/range_)**3)
            else:
                return nugget + sill

    # 构建方程矩阵
    A = np.zeros((n+1, n+1))
    for i in range(n):
        for j in range(n):
            A[i, j] = variogram(distance_matrix[i, j], model, nugget, sill, range_)
        A[i, n] = 1
        A[n, i] = 1

    # 构建值向量
    b = np.zeros(n+1)
    for i in range(n):
        b[i] = variogram(np.linalg.norm(points[i] - target_point), model, nugget, sill, range_)
    b[n] = 1

    # 解方程
    weights = np.linalg.solve(A, b)

    # 计算插值结果
    interpolated_value = np.dot(values, weights[:-1])
    return interpolated_value


if __name__ == '__main__':
    pass
    # Data = pd.read_excel("Interpolate.xlsx")
    # Points = Data.loc[:, ['经度','纬度']].values
    # Values = Data.loc[:, ['值']].values
    #
    # # 普通克里金（球面函数模型）插值
    # KD = gma.smc.Interpolate.Kriging(Points, Values, Resolution = 0.05,
    #                              VariogramModel = 'Spherical',
    #                              VariogramParameters = None,
    #                              KMethod = 'Ordinary',
    #                              InProjection = 'EPSG:4326')
    #
    # # 泛克里金类似，这里不做示例
    #
    # gma.rasp.WriteRaster(r'.\gma_OKriging.tif',
    #                  KD.Data,
    #                  Projection = 'WGS84',
    #                  Transform = KD.Transform,
    #                  DataType = 'Float32')